Revised Interoceptive Accuracy Scale (IAS-R)

Analysis

Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Study 1 - Reanalysis of Murphy et al. (2020)

Reanalyze their data to confirm the factor structure of the IAS.

Data Preparation

Code
# https://osf.io/3m5nh
df1 <- haven::read_sav("../data/Murphy2020/Study 1.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other")))) 

# https://osf.io/3m5nh
df2 <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other")))) 

# Campos (2022) - Study 1 (https://osf.io/j6ef3)
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |> 
  rename(Heart=IAS1, Hungry=IAS2, Breathing=IAS3, Thirsty=IAS4,
         Urinate=IAS5, Defecate=IAS6, Taste=IAS7, Vomit=IAS8, Sneeze=IAS9,
         Cough=IAS10, Temp=IAS11, Sex_arousal=IAS12, Wind=IAS13, Burp=IAS14,
         Muscles=IAS15, Bruise=IAS16, Pain=IAS17, Blood_Sugar=IAS18,
         Affective_touch=IAS19, Tickle=IAS20, Itch=IAS21)

# https://osf.io/3eztd
df4 <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |> 
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |> 
  rename(Age=age, 
         Heart=Heart, Hungry=HUNGR, Breathing=BREAT, Thirsty=Thirs,
         Urinate=URINA, Defecate=Defec, Taste=TASTE, Vomit=VOMIT, Sneeze=Sneez,
         Cough=COUGH, Temp=TEMPE, Sex_arousal=SEXAR, Wind=WIND, Burp=Burp,
         Muscles=MUSCL, Bruise=Bruis, Pain=PAIN, Blood_Sugar=BloSu,
         Affective_touch=Touch, Itch=ITCH)  # No tickle because same Chinese character

# https://osf.io/3eztd
df5 <- haven::read_sav("../data/Lin2023/Study 2.sav") |> 
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |> 
  rename(Heart=Heart, Hungry=HUNGR, Breathing=BREAT, Thirsty=Thirs,
         Urinate=URINA, Defecate=Defec, Taste=TASTE, Vomit=VOMIT, Sneeze=Sneez,
         Cough=COUGH, Temp=TEMPE, Sex_arousal=SEXAR, Wind=WIND, Burp=Burp,
         Muscles=MUSCL, Bruise=Bruis, Pain=PAIN, Blood_Sugar=BloSu,
         Affective_touch=Touch, Itch=ITCH)  # No tickle because same Chinese character

df6 <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(Heart=IAS_1, Hungry=IAS_2, Breathing=IAS_3, Thirsty=IAS_4, 
         Urinate=IAS_5, Defecate=IAS_6, Taste=IAS_7, Vomit=IAS_8, Sneeze=IAS_9, 
         Cough=IAS_10, Temp=IAS_11, Sex_arousal=IAS_12, Wind=IAS_13, Burp=IAS_14,
         Muscles=IAS_15, Bruise=IAS_16, Pain=IAS_17, Blood_Sugar=IAS_18,
         Affective_touch=IAS_19, Tickle=IAS_20, Itch=IAS_21)

df7 <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(Gender=Sex, Heart=Item_IAS_Interoception_1, Hungry=Item_IAS_Interoception_2, 
         Breathing=Item_IAS_Interoception_3, Thirsty=Item_IAS_Interoception_4, 
         Temp=Item_IAS_Interoception_5, Sex_arousal=Item_IAS_Interoception_6,
         Urinate=Item_IAS_Elimination_1, Defecate=Item_IAS_Elimination_2,
         Vomit=Item_IAS_Elimination_3, Wind=Item_IAS_Expulsion_1, 
         Burp=Item_IAS_Expulsion_2, Sneeze=Item_IAS_Expulsion_3,
         Muscles=Item_IAS_Nociception_1, Bruise=Item_IAS_Nociception_2,
         Pain=Item_IAS_Nociception_3, Affective_touch=Item_IAS_Skin_1, 
         Tickle=Item_IAS_Skin_2, Itch=Item_IAS_Skin_3) |>
  filter(!is.na(Urinate))

Participants

  • Sample 1: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
  • Sample 2: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
  • Sample 3: Data from Campos’ (2022) study 1, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
  • Sample 4: Data from Lin’s (2023) study 1 and 3, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
  • Sample 5: Data from Lin’s (2023) study 2, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
  • Sample 6: New dataset from Makowski (under preparation), included 251 participants (Mean age = 37.5, SD = 13.2, range: [17, 76]; Gender: 53.0% women, 42.6% men, 2.79% non-binary, 1.59% missing; Education: Bachelor, 39.84%; Doctorate, 4.78%; High School, 31.47%; Master, 19.12%; missing, 0.40%; Other, 4.38%).
  • Sample 7: New dataset from Makowski (under preparation), included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).

Total N = 3743.

Descriptive

Distribution

Code
vars <- c("Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
          "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch")
  
dens1 <- estimate_density(select(df1, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 1")
dens2 <- estimate_density(select(df2, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(setdiff(vars, "Tickle"))), method="kernSmooth") |> 
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(setdiff(vars, "Tickle"))), method="kernSmooth") |> 
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 6")
dens7 <- estimate_density(select(df7, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method="kernSmooth") |> 
  mutate(Sample = "Sample 7")

rbind(dens1, dens2, dens3, dens4, dens5, dens6, dens7) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  theme_minimal() +
  labs(title = "Distribution of Responses", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free", nrow=4)

Code
data1 <- normalize(select(df1, all_of(dens1$Parameter)), range=c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range=c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range=c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range=c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range=c(1, 5))
data6 <- select(df6, all_of(dens6$Parameter))
data7 <- select(df7, all_of(dens7$Parameter))

data_all <- rbind(data1, data2, data3, 
                  mutate(data4, Tickle=NA), mutate(data5, Tickle=NA), 
                  data6, mutate(data7, Taste=NA, Cough=NA, Blood_Sugar=NA)) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

Code
make_correlation <- function(df) {
  correlation::correlation(df, redundant=TRUE) |> 
    correlation::cor_sort() |> 
    correlation::cor_lower() |> 
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only=TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x=Parameter1, y=Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks=FALSE)) +
    theme_minimal() +
    theme(legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1))
}

make_correlation(data_all)

Code
make_correlation(data1)

Code
make_correlation(data2)

Code
make_correlation(data3)

Code
make_correlation(data4)

Code
make_correlation(data5)

Code
make_correlation(data6)

Code
make_correlation(data7)

EGA

See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

Code
uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.306

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.222
 Urinate Defecate 0.219
  Sneeze    Cough 0.204
Code
uva$keep_remove
$keep
[1] "Itch"

$remove
[1] "Tickle"
Code
uva <- EGAnet::UVA(data = data1, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"
Code
uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208
Code
uva$keep_remove
$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"
Code
uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data7, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209
Code
uva$keep_remove
NULL

Item Stability

Code
plots <- list()
for(model in c("glasso", "TMFG")) {
  for(algo in c("walktrap", "louvain")) {
    for(type in c("ega", "ega.fit")) {  # "hierega"
      print(type)
      ega <- EGAnet::bootEGA(
        data = data_all,
        seed=123,
        model=model,
        algorithm=algo,
        EGA.type=type,
        type="resampling",
        plot.itemStability=FALSE,
        verbose=FALSE)
      plots[[length(plots) + 1]] <- plot(ega) + 
        labs(title = paste0(model, " (", algo, ") - ", type))
      }
    }
}
[1] "ega"
[1] "ega.fit"
[1] "ega"
[1] "ega.fit"
[1] "ega"
[1] "ega.fit"
[1] "ega"
[1] "ega.fit"
Code
length(plots)
[1] 8
Code
patchwork::wrap_plots(plots, nrow = 4)

Exploratory Factor Analysis (EFA)

How Many Factors

Code
n <- parameters::n_factors(data_all, n_max = 12)

plot(n) 

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4, 5 ), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
Code
n <- parameters::n_factors(data1, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data2, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data3, n_max = 12)

plot(n) 

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
4 BIC BIC
Code
n <- parameters::n_factors(data4, n_max = 12)

plot(n) 

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data5, n_max = 12)

plot(n) 

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 BIC BIC
Code
n <- parameters::n_factors(data6, n_max = 12)

plot(n) 

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data7, n_max = 12)

plot(n) 

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 BIC BIC
4 BIC (adjusted) BIC

EFA Models

Code
efa5 <- parameters::factor_analysis(data_all, n=5, rotation = "oblimin", sort = TRUE)
Loading required namespace: GPArotation
Code
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR3 MR2 MR1 MR5 MR4 Complexity Uniqueness
Cough 0.72 -5.89e-04 -0.06 -5.43e-03 0.03 1.02 0.50
Sneeze 0.67 -0.03 0.04 0.06 -0.05 1.04 0.53
Burp 0.55 0.13 0.06 -0.03 -0.02 1.15 0.59
Wind 0.48 0.11 -9.06e-03 -1.90e-03 0.06 1.15 0.68
Vomit 0.35 0.03 0.24 0.14 -0.09 2.32 0.68
Itch 0.05 0.71 -0.10 -0.06 0.12 1.12 0.46
Tickle 0.08 0.68 -0.03 0.06 4.30e-03 1.05 0.45
Bruise 8.13e-04 0.45 0.22 0.06 -0.08 1.56 0.68
Blood_Sugar 0.06 0.36 0.14 0.04 -0.06 1.43 0.79
Affective_touch -0.05 0.36 0.34 0.18 -0.16 2.91 0.64
Muscles 0.03 0.25 0.21 0.17 0.07 2.95 0.72
Taste 0.16 0.18 0.17 0.13 0.05 3.96 0.75
Urinate 0.10 -0.08 0.55 0.01 0.20 1.39 0.55
Defecate 0.22 -0.04 0.46 -0.01 0.15 1.69 0.60
Pain -0.01 0.22 0.28 0.11 0.23 3.22 0.65
Heart 1.28e-03 -0.03 -0.01 0.66 -0.01 1.01 0.59
Breathing 0.07 7.82e-03 -0.08 0.56 0.14 1.20 0.61
Sex_arousal 0.10 0.07 0.18 0.25 0.12 2.87 0.73
Thirsty -1.95e-04 0.05 0.11 6.26e-03 0.62 1.07 0.53
Hungry 0.01 0.04 -5.98e-03 0.12 0.60 1.08 0.55
Temp 0.11 0.12 0.14 0.15 0.27 3.12 0.68

The 5 latent factors (oblimin rotation) accounted for 38.36% of the total variance of the original data (MR3 = 10.15%, MR2 = 9.46%, MR1 = 6.49%, MR5 = 6.20%, MR4 = 6.05%).

Code
efa4 <- parameters::factor_analysis(data2, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR3 MR4 Complexity Uniqueness
Urinate 0.69 -0.18 0.13 -0.04 1.22 0.54
Defecate 0.68 -0.13 0.04 0.06 1.10 0.55
Breathing 0.56 0.03 0.06 7.96e-03 1.03 0.63
Hungry 0.55 0.15 -7.15e-03 0.02 1.16 0.60
Thirsty 0.52 0.05 0.14 -0.01 1.16 0.62
Sex_arousal 0.50 0.14 -0.04 0.20 1.49 0.56
Temp 0.46 0.23 0.16 -0.01 1.77 0.52
Pain 0.46 0.43 -0.05 -0.03 2.03 0.49
Taste 0.40 0.21 -0.07 0.12 1.77 0.69
Heart 0.40 0.06 -0.04 0.15 1.36 0.76
Muscles 0.35 0.24 0.07 0.09 2.07 0.65
Vomit 0.32 0.06 0.23 0.11 2.20 0.68
Tickle 0.02 0.69 0.01 0.04 1.01 0.47
Itch -0.02 0.67 0.18 0.03 1.16 0.41
Bruise 0.07 0.44 0.15 0.06 1.32 0.65
Blood_Sugar -0.08 0.38 0.16 0.09 1.61 0.77
Affective_touch 0.38 0.38 -0.13 0.05 2.26 0.63
Sneeze 0.05 0.05 0.71 0.01 1.02 0.41
Cough 0.04 0.02 0.71 0.04 1.01 0.43
Wind 7.52e-03 -0.04 -5.18e-03 0.94 1.00 0.15
Burp -0.01 0.27 0.17 0.47 1.91 0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Code
efa4 <- parameters::factor_analysis(data2, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR3 MR4 Complexity Uniqueness
Urinate 0.69 -0.18 0.13 -0.04 1.22 0.54
Defecate 0.68 -0.13 0.04 0.06 1.10 0.55
Breathing 0.56 0.03 0.06 7.96e-03 1.03 0.63
Hungry 0.55 0.15 -7.15e-03 0.02 1.16 0.60
Thirsty 0.52 0.05 0.14 -0.01 1.16 0.62
Sex_arousal 0.50 0.14 -0.04 0.20 1.49 0.56
Temp 0.46 0.23 0.16 -0.01 1.77 0.52
Pain 0.46 0.43 -0.05 -0.03 2.03 0.49
Taste 0.40 0.21 -0.07 0.12 1.77 0.69
Heart 0.40 0.06 -0.04 0.15 1.36 0.76
Muscles 0.35 0.24 0.07 0.09 2.07 0.65
Vomit 0.32 0.06 0.23 0.11 2.20 0.68
Tickle 0.02 0.69 0.01 0.04 1.01 0.47
Itch -0.02 0.67 0.18 0.03 1.16 0.41
Bruise 0.07 0.44 0.15 0.06 1.32 0.65
Blood_Sugar -0.08 0.38 0.16 0.09 1.61 0.77
Affective_touch 0.38 0.38 -0.13 0.05 2.26 0.63
Sneeze 0.05 0.05 0.71 0.01 1.02 0.41
Cough 0.04 0.02 0.71 0.04 1.01 0.43
Wind 7.52e-03 -0.04 -5.18e-03 0.94 1.00 0.15
Burp -0.01 0.27 0.17 0.47 1.91 0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Code
efa4 <- parameters::factor_analysis(data3, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Cough 0.88 -0.08 -0.04 0.11 1.05 0.29
Sneeze 0.74 0.05 1.42e-03 -0.03 1.01 0.41
Vomit 0.55 0.12 0.04 -0.09 1.17 0.58
Burp 0.55 0.08 0.22 -0.13 1.47 0.48
Wind 0.53 0.10 0.21 -0.25 1.84 0.48
Temp 0.46 0.34 -0.06 0.09 1.97 0.52
Taste 0.29 0.18 0.10 0.12 2.37 0.74
Hungry -8.85e-03 0.70 -0.01 -0.09 1.03 0.54
Urinate 0.08 0.66 -0.12 0.11 1.15 0.52
Thirsty -0.08 0.64 -0.03 0.14 1.14 0.62
Defecate 0.10 0.63 -3.50e-05 -0.07 1.08 0.53
Breathing -0.10 0.49 0.28 -0.07 1.76 0.66
Sex_arousal 0.17 0.44 0.17 -0.22 2.20 0.58
Pain 0.21 0.41 0.08 0.18 2.00 0.58
Heart 0.05 0.28 0.24 0.08 2.23 0.77
Muscles 0.24 0.27 0.20 0.19 3.66 0.61
Tickle 0.02 -0.05 0.84 -0.04 1.01 0.31
Itch 0.03 0.02 0.74 0.11 1.05 0.39
Affective_touch 0.08 0.23 0.36 0.12 2.08 0.67
Bruise 0.17 0.07 0.29 0.45 2.08 0.55
Blood_Sugar 4.11e-03 0.11 0.35 0.43 2.07 0.60

The 4 latent factors (oblimin rotation) accounted for 45.69% of the total variance of the original data (MR1 = 15.50%, MR3 = 15.46%, MR2 = 11.02%, MR4 = 3.71%).

Code
efa5 <- parameters::factor_analysis(data4, n=5, rotation = "oblimin", sort = TRUE)
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR5 MR4 MR2 Complexity Uniqueness
Cough 0.67 -0.07 8.90e-03 0.06 0.03 1.04 0.54
Sneeze 0.63 -0.04 0.11 0.05 -0.08 1.11 0.54
Burp 0.58 0.18 -0.03 -0.09 0.03 1.25 0.59
Wind 0.54 0.08 -0.06 -0.03 0.11 1.17 0.67
Vomit 0.33 0.06 0.21 0.15 -0.19 3.02 0.70
Taste 0.24 0.19 -4.60e-03 0.06 0.18 3.00 0.77
Bruise -0.03 0.63 0.05 0.05 -0.03 1.03 0.58
Affective_touch 0.07 0.53 0.03 -0.01 -0.04 1.05 0.67
Blood_Sugar 0.07 0.43 0.15 0.10 -0.18 1.82 0.67
Itch 0.11 0.41 -0.01 0.05 0.15 1.45 0.71
Muscles 0.07 0.40 -0.05 0.02 0.19 1.56 0.75
Urinate -0.04 0.06 0.64 6.75e-03 0.06 1.04 0.55
Defecate 0.11 0.02 0.59 -0.03 0.02 1.08 0.57
Heart -0.04 0.04 -0.03 0.66 -0.04 1.02 0.59
Breathing 0.09 -2.40e-03 -0.02 0.51 0.12 1.17 0.65
Sex_arousal 0.07 0.18 0.07 0.24 0.13 2.86 0.78
Thirsty 0.05 -0.05 0.24 0.03 0.49 1.51 0.60
Hungry 0.09 -0.10 0.15 0.21 0.38 2.25 0.67
Temp 0.09 0.17 7.51e-03 0.11 0.36 1.81 0.72
Pain -3.56e-03 0.30 0.14 0.04 0.33 2.38 0.67

The 5 latent factors (oblimin rotation) accounted for 35.08% of the total variance of the original data (MR1 = 10.08%, MR3 = 8.61%, MR5 = 5.94%, MR4 = 5.33%, MR2 = 5.12%).

Code
efa5 <- parameters::factor_analysis(data5, n=5, rotation = "oblimin", sort = TRUE)
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 MR5 Complexity Uniqueness
Burp 0.69 0.06 0.03 -0.10 0.03 1.06 0.50
Wind 0.62 0.01 0.03 -0.02 -0.08 1.04 0.64
Cough 0.54 -0.06 0.01 0.17 0.11 1.31 0.58
Sneeze 0.45 0.02 -0.07 0.21 0.17 1.78 0.61
Bruise -0.04 0.60 -0.01 0.02 0.06 1.03 0.63
Affective_touch 0.08 0.53 -0.04 -0.08 0.06 1.12 0.67
Blood_Sugar 0.01 0.48 -0.09 0.08 0.19 1.45 0.68
Itch 0.18 0.38 0.08 0.12 -0.02 1.77 0.69
Muscles 0.09 0.38 0.19 3.00e-03 -0.07 1.68 0.77
Thirsty 0.03 -0.17 0.57 -0.01 0.22 1.50 0.55
Temp 2.56e-03 0.18 0.49 0.03 -0.08 1.34 0.71
Hungry 0.02 -0.02 0.48 0.16 0.09 1.32 0.65
Taste 0.19 0.22 0.32 0.05 -0.12 2.84 0.73
Pain 0.10 0.24 0.28 2.96e-03 0.03 2.26 0.77
Breathing 0.04 -0.03 0.03 0.67 -0.07 1.03 0.54
Heart -0.08 0.05 5.68e-03 0.58 0.06 1.08 0.66
Sex_arousal 0.11 0.14 0.23 0.26 -0.06 3.08 0.76
Urinate -0.04 0.11 0.11 -0.01 0.61 1.14 0.54
Defecate 0.14 0.03 0.03 -0.02 0.60 1.11 0.53
Vomit 0.15 0.15 -0.04 0.22 0.23 3.56 0.74

The 5 latent factors (oblimin rotation) accounted for 35.25% of the total variance of the original data (MR1 = 8.99%, MR4 = 8.07%, MR2 = 6.47%, MR3 = 5.97%, MR5 = 5.75%).

Code
efa4 <- parameters::factor_analysis(data5, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 Complexity Uniqueness
Burp 0.65 0.07 6.06e-03 -0.10 1.07 0.56
Cough 0.63 -0.06 0.03 0.11 1.09 0.56
Sneeze 0.57 0.02 -8.72e-03 0.13 1.10 0.61
Wind 0.54 0.02 -0.03 1.35e-04 1.01 0.71
Defecate 0.34 0.11 0.28 -0.13 2.48 0.70
Vomit 0.27 0.18 0.06 0.13 2.39 0.77
Bruise -0.06 0.64 -0.04 0.01 1.02 0.63
Affective_touch 0.08 0.55 -0.07 -0.08 1.11 0.67
Blood_Sugar 0.09 0.50 -0.03 0.02 1.07 0.71
Muscles -2.84e-04 0.41 0.12 0.05 1.19 0.79
Itch 0.15 0.40 0.02 0.13 1.52 0.69
Pain 0.05 0.29 0.26 0.02 2.07 0.78
Taste 0.09 0.24 0.21 0.11 2.76 0.79
Thirsty 0.03 -0.10 0.71 -0.04 1.05 0.51
Hungry -4.08e-03 0.04 0.51 0.17 1.23 0.65
Temp -0.09 0.22 0.39 0.10 1.87 0.77
Urinate 0.19 0.18 0.33 -0.11 2.49 0.73
Breathing 0.06 -0.03 -0.01 0.70 1.02 0.50
Heart -0.01 0.07 0.04 0.50 1.06 0.72
Sex_arousal 0.06 0.17 0.17 0.30 2.35 0.76

The 4 latent factors (oblimin rotation) accounted for 32.02% of the total variance of the original data (MR1 = 9.96%, MR4 = 9.20%, MR2 = 7.36%, MR3 = 5.49%).

Model Construction

Structure

Code
m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"


m7 <- "
# ----
Sustenance =~ Hungry + Thirsty
Elimination =~ Urinate + Defecate
Skin =~ Itch + Bruise
Nociception =~ Muscles + Pain
Affect =~ Breathing + Heart
# ----
Expulsion =~ Cough + Sneeze + Wind + Burp

# --- 
# Ambiguous: Temp + Vomit + Affective_touch + Sex_arousal + Taste
# Discard: Tickle + Blood_Sugar
" 

m8 <- "
# ----
Sustenance =~ Hungry + Thirsty
Elimination =~ Urinate + Defecate
Skin =~ Itch + Bruise
Nociception =~ Muscles + Pain
Affect =~ Breathing + Heart
# ----
Expulsion_Sudden =~ Cough + Sneeze 
Expulsion_Gastric =~ Wind + Burp
" 
# Ambiguous: Temp + Vomit + Affective_touch + Sex_arousal + Taste
# Discard: Tickle + Blood_Sugar

m1_all <- lavaan::cfa(m1, data = data_all)
m7_all <- lavaan::cfa(m7, data = data_all)
m8_all <- lavaan::cfa(m8, data = data_all)

anova(m1_all, m7_all, m8_all) |> 
  parameters::parameters() |> 
  display()
Parameter AIC BIC Chi2 df p
m8_all -14218.06 -13919.70 239.58 56
m7_all -14036.69 -13774.87 432.95 62 < .001
m1_all -12372.54 -12202.05 2127.10 77 < .001

Anova Table (Type 1 tests)

Higher Order Factors

Code
m8h1 <- paste0(m8, "
              Interoception =~ Sustenance + Elimination + Skin + Nociception + Affect + Expulsion_Sudden + Expulsion_Gastric
              ")

m8h3 <- paste0(m8, "
              Valenced =~ Nociception + Affect + Skin
              Expulsion =~ Expulsion_Sudden + Expulsion_Gastric 
              Homeostasis =~ Sustenance + Elimination
              ")

m8h1_all <- lavaan::cfa(m8h1, data = data_all)
m8h3_all <- lavaan::cfa(m8h3, data = data_all)
anova(m8_all, m8h1_all, m8h3_all) |> 
  parameters::parameters() |> 
  display()
Parameter AIC BIC Chi2 df p
m8_all -14218.06 -13919.70 239.58 56
m8h3_all -14104.40 -13873.02 375.24 67 < .001
m8h1_all -13903.22 -13690.11 582.42 70 < .001

Anova Table (Type 1 tests)

Model Performance

Code
rbind(
  performance::performance(m8_all, metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |> 
    mutate(Sample = "All"),
  performance::performance(update(m8_all, data=data1), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |> 
    mutate(Sample = "Sample 1"),
  performance::performance(update(m8_all, data=data2), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m8_all, data=data3), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m8_all, data=data4), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m8_all, data=data5), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m8_all, data=data6), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(lavaan::cfa(str_remove(m8, fixed("\nExpulsion_Sudden =~ Cough + Sneeze")), data=data7), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7")
) |> 
  display()
Chi2 RMSEA CFI SRMR Sample
239.58 0.03 0.98 0.02 All
88.22 0.04 0.98 0.03 Sample 1
66.50 0.02 0.99 0.03 Sample 2
114.03 0.04 0.98 0.03 Sample 3
115.89 0.03 0.98 0.02 Sample 4
83.61 0.03 0.98 0.03 Sample 5
91.47 0.05 0.97 0.04 Sample 6
82.59 0.05 0.97 0.03 Sample 7

No evidence for higher order factors. The unidimensionality might be a byproduct of the similarities in phrasing of the items.

Study 2 - Correlates

Code
df1 <- cbind(df1, as.data.frame(predict(lavaan::cfa(m8, data = data1)))) |> 
  mutate(Sample = "Sample 1")
df2 <- cbind(df2, as.data.frame(predict(lavaan::cfa(m8, data = data2)))) |> 
  mutate(Sample = "Sample 2")
df3 <- cbind(df3, as.data.frame(predict(lavaan::cfa(m8, data = data3)))) |> 
  mutate(Sample = "Sample 3")
df4 <- cbind(df4, as.data.frame(predict(lavaan::cfa(m8, data = data4)))) |> 
  mutate(Sample = "Sample 4")
df5 <- cbind(df5, as.data.frame(predict(lavaan::cfa(m8, data = data5)))) |> 
  mutate(Sample = "Sample 5")
df6 <- cbind(df6, as.data.frame(predict(lavaan::cfa(m8, data = data6)))) |> 
  mutate(Sample = "Sample 6")
df7 <- cbind(df7, as.data.frame(predict(lavaan::cfa(str_remove(m8, fixed("\nExpulsion_Sudden =~ Cough + Sneeze")), data = data7)))) |> 
  mutate(Sample = "Sample 7", Expulsion_Sudden = NA)

vars <- names(as.data.frame(predict(m8_all)))

Demographics

Code
df_demo <- rbind(
  select(df1, Sample, Age, Gender, all_of(vars)),
  select(df2, Sample, Age, Gender, all_of(vars)),
  select(df3, Sample, Age, Gender, all_of(vars)),
  select(df4, Sample, Age, Gender, all_of(vars)),
  select(df5, Sample, Age, Gender, all_of(vars)),
  select(df6, Sample, Age, Gender, all_of(vars)),
  select(df7, Sample, Age, Gender, all_of(vars))
) 

make_lm_data <- function(df, vars, predictor="Age") {
  dat <- data.frame()
  for(resp in vars) {
      if(length(unique(df$Sample)) == 1) {
        m <- lm(as.formula(paste0(resp, " ~ ", predictor)), data = df)
        param <- parameters::parameters(m)[2, ] |> 
          mutate(Sample = unique(df$Sample)) |> 
          select(-Parameter)
      } else {
        m <- lm(as.formula(paste0(resp, " ~ Sample / ", predictor)), data = df)
        param <- parameters::parameters(m) |> 
          filter(str_detect(Parameter, ":")) |> 
          separate(Parameter, c("Predictor", "Sample"), sep = "Sample ") |> 
          mutate(Sample = paste0("Sample ", str_remove(Sample, fixed(paste0(":", predictor)))))
      }
    
      dat <- param |> 
        mutate(Predictor = predictor,
               Repsonse = resp,
               t_low = CI_low / SE,
               t_high = CI_high / SE) |> 
        rbind(dat)
  }
  dat
}

data <- rbind(
  make_lm_data(df_demo, vars, predictor="Age"),
  make_lm_data(filter(df_demo, Gender %in% c("Male", "Female")), vars, predictor="Gender")
)

BPQ

Code
x <- intersect(names(df3)[str_detect(names(df3), "BPQ")], 
               names(df4)[str_starts(names(df4), "BPQ")])

df_bpq <- rbind(
  select(df3, Sample, all_of(x), all_of(vars)),
  select(df4, Sample, all_of(x), all_of(vars))
) 

# estimate_density(select(df_bpq, Sample, starts_with("BPQ")), at="Sample")  |> 
#   plot() +
#   facet_grid(~Sample)

plot(parameters::n_components(select(df_bpq, starts_with("BPQ"))))

Code
df_bpq$BPQ_PCA1 <- predict(parameters::principal_components(select(df_bpq, starts_with("BPQ")), n=1))$PC1

data <- rbind(data,
              make_lm_data(df_bpq, vars, predictor="BPQ_PCA1"))

# Why opposite effects?

TAS

Code
x <- intersect(names(df3)[str_detect(names(df3), "TAS")], 
               names(df5)[str_starts(names(df5), "TAS")])

PHQ

Code
names(df5)[str_detect(names(df5), "PHQ")]
 [1] "PHQ15_sum" "PHQ9_sum"  "PHQ9_1"    "PHQ9_2"    "PHQ9_3"    "PHQ9_4"   
 [7] "PHQ9_5"    "PHQ9_6"    "PHQ9_7"    "PHQ9_8"    "PHQ15_1"   "PHQ15_2"  
[13] "PHQ15_3"   "PHQ15_4"   "PHQ15_5"   "PHQ15_6"   "PHQ15_7"   "PHQ15_8"  
[19] "PHQ15_9"   "PHQ15_10"  "PHQ15_11"  "PHQ15_12"  "PHQ15_13" 
Code
data <- rbind(data,
              make_lm_data(df5, vars, predictor="PHQ9_sum"),
              make_lm_data(df5, vars, predictor="PHQ15_sum"))

BDI

Code
# names(df1)
# names(df2)
# names(df3)  # TAS, BPQ
# names(df4)  # BPQ
# names(df5)  # PHQ, TAS
# names(df6)  # BDI, STAI
# names(df7)  # MAIA

data <- rbind(data,
              make_lm_data(df6, vars, predictor="BDI2_Total"))

Summary

Code
data |> 
  mutate(sig = ifelse(p < .01, "p < .01", ifelse(p < .001, "p < .001", "N.S."))) |> 
  ggplot(aes(x=Predictor)) +
  geom_hline(yintercept=0, linetype="dashed") +
  geom_pointrange(aes(y=t, ymin=t_low, ymax=t_high, color=sig), position=position_dodge2(width=0.3)) +
  facet_grid(~Repsonse, scales = "free_y") +
  coord_flip()

Discussion

Benefits of the IAS: - Straightforward and sensation-centered items

Recommendations: - Remove Itch (redundant + issue in Chinese) - Use analog scales

Limitations: - Not much clear theorethical or empirical structure (small grouping of items) - Limited variability (clear mode at 4/5) - Ambiguous items which grouping depends on the context (and its awareness) - E.g., heart beating fast, vomit when nauseaous - Few items for some modalities (e.g., 1 for heart) - Positive phrasing: benefits but also might exacerbate positivity bias (and thus unidimensionality)

Need for context-specific items (cross-modal when possible, i.E., cardioception, respiroception, etc.).

New Scale: Multimodal Interoceptive Sensitivity Scale (MISS)